%pylab inline
import matplotlib.pyplot as plt

from ipywidgets import *

Populating the interactive namespace from numpy and matplotlib

Kausztika. Miért csillog a Balaton vízfelszíne?

Kiegészítés: Cserti József: Kausztikák a fizikában

Lásd még: 2018. évi Ortvay verseny 18. feladata

def reflection(be,nn):
    '''
    be: a bejovo fenysugar iranyvektora (normalt)
    nn: a felulet normalvektora (normalt)
    r: a kimeno fenysugar iranyvektora (normalt)
    
    '''
    #be=be/sqrt((sum(be**2)))  # normalizing

    ki=be-2*nn*dot(be,nn)
    return(ki)
    

A felület alakja:

$y(x)=A*\sin(2\pi x/\lambda)$

A Balaton vízfelszínének egy modellje

def ref_rays(be,x0,param):
    '''
    a felulet fv: y(x)=A*sin(2\pi x/\lambda)
    
    '''
    A=1
    lam=20
    A,lam = param
    
    y0=A*sin(2*pi*x0/lam)
    der=A*2*pi/lam*cos(2*pi*x0/lam)
    nn=array([der,-1])    
    normal=nn/sqrt((sum(nn**2))) # normalizing
    
    ki=reflection(be,normal)
    
    return(y0,ki)
    

def gorbe(p,A,k,theta):
    '''
    a kausztika parameteres egyenlete, 'p' a parameter
    'theta' a fuggolegesen lefele mutato egysegvektor iranyanak szoge a fuggolegeshez kepest.
    
    '''
    
    xx1=(2 + A**2 * k**2 + A**2 *k**2 *cos(2*k* p))/sin(k*p)*(cos(theta)+A*k*cos(k* p)* sin(theta))
    xx2=2*A*k*cos(k* p)*cos(theta)-sin(theta)+A**2 * k**2 *(cos(k* p))**2 *sin(theta)
    nev=4*k**2 *(A+A**3 * k**2 *(cos(k* p))**2)
    xk=p + xx1*xx2/nev
    
    yy1=A*k*cos(k* p)*sin(theta)+cos(theta)
    yy2=(-1+A**2 * k**2 *(cos(k* p))**2)* cos(theta) -2* A*k*cos(k* p)* sin(theta)
    yk=A* sin(k* p) + (2+A**2 * k**2 + A**2 *k**2 *cos(2*k* p))/sin(k*p)*yy1*yy2/nev
    
    
    return(xk,yk)

def rajz_Balaton_kausztika_eq(A,k,theta):
    
    #A,k=(1,0.5)
    #alpha = 65*pi/180, itt  alpha= pi/2 -theta, theta a beesesi szog a fuggolegeshez merve

    Np=5000
    pmax=50
    pmin=-pmax
    p=linspace(pmin,pmax,Np)
    yfelszin=A*sin(k*p)
    p=linspace(pmin,pmax,Np)
    xk,yk=gorbe(p,A,k,theta)

    #figsize(12,12)
    #ax=subplot(aspect='equal')
    #ax=subplot(111)

    #plot(p,yfelszin, color='b')
    plot(xk,yk,'b.');

    #xlim(-15,15)
    #ylim(-15,15);
    #ax.set_axis_off();
    grid();
    

theta=10*pi/180
be=array([sin(theta),-cos(theta)])
#be=be/sqrt((sum(be**2)))  # normalizing

A,lam=(1,20)
param= (A,lam)

Np,xmin,xmax=(100,-15,45)
ymax=20*A       
    
xlist=linspace(xmin,xmax,Np)
ylist=A*sin(2*pi*xlist/lam)

figsize(12,8)
#ax=subplot(aspect='equal')
ax=subplot(111)


plot(xlist,ylist,color='k',lw=3)

benyil=0.25*ymax*be
arrow(lam/3,ymax,benyil[0],benyil[1], head_width=0.5, head_length=0.5, fc='b', ec='g',lw=7)

for x0 in xlist:
    y0,rr=ref_rays(be,x0,param)

    y1=ymax
    x1=x0+rr[0]/rr[1]*(y1-y0)

    plot([x0,x1],[y0,y1],color='r');
        
k=2*pi/lam
#theta = 10*pi/180
    
rajz_Balaton_kausztika_eq(A,k,theta)


xlim(xmin,xmax)
ylim(-2.5*A,ymax)
title('Visszavert fénysugarak (piros), bejövő (zöld)\n', fontsize=18);
#ax.set_axis_off();

def rajz_Balaton(theta,A,lam,Np,xmeret,ymax):
    
    theta=theta*pi/180
    be=array([sin(theta),-cos(theta)])
    #be=be/sqrt((sum(be**2)))  # normalizing

    #A,lam=(1,20)
    param= (A,lam)

    xmin,xmax = xmeret
    #xmin,xmax=(-15,45)
    #ymax=10*A
    
    xlist=linspace(xmin,xmax,Np)
    ylist=A*sin(2*pi*xlist/lam)

    figsize(10,6)
    #ax=subplot(aspect='equal')
    ax=subplot(111)

    plot(xlist,ylist,color='k',lw=3)

    benyil=0.25*ymax*be
    arrow(0,ymax,benyil[0],benyil[1], head_width=0.5, head_length=0.5, fc='b', ec='b',lw=7)

    for x0 in xlist:
        y0,rr=ref_rays(be,x0,param)

        y1=ymax
        x1=x0+rr[0]/rr[1]*(y1-y0)

        plot([x0,x1],[y0,y1],color='r');
        
    xlim(xmin,xmax)
    ylim(-1.5*A,ymax)
    title('Visszavert fénysugarak (piros), bejövő (kék)\n', fontsize=18);

    rajz_Balaton_kausztika_eq(A,k,theta)

    ax.set_axis_off();

xmeret=(-35,25)  # xmin,xmax,ymax
A,lam=(1,20)
              
@interact(theta=(-90,90,5),Np=(50,300,10),ymax=(10,100,10))
def play(theta=0,Np=100,ymax=25):
    
    rajz_Balaton(theta,A,lam,Np,xmeret,ymax)
    show();